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Abstract: In this paper the multi-dimensional random walk models governed 
by distributed fractional order differential equations and multi-term fractional 
order differential equations are constructed. The scaling limits of these random 
walks to a diffusion process in the sense of distributions is proved. 



1. Introduction 



In this paper we construct new random walks connected with fractional order differ- 
ential equations. Namely, the governing equations corresponding to the constructed 
random walks are multi-term or distributed fractional order differential equations. 
Nowadays the connection between random walk and fractional order dynamics is 
well known, see, for instance [H, Eli Ell 0, 38 1. A number of constructive random 



walk models governed by fractional differential equations in the one-dimensional 
case were studied by Gillis, et al. [T3|, Chechkin, et al. 0, Gorenflo, et al. [3 El, 



and in the n-dimensional case by Umarov [35(, Umarov, et al. [361 ] . Andries, et al. 

The governing equation in these studies depends on parameters /3 G (0, 1] and 
a G (0, 2], and is given by the fractional order differential equation 



(1) V l3 u(t,x)=D%u(t,x),t>0,x<E 



where X" 9 is the time-fractional derivative in some sense, and Dq , < a < 2, is the 
pseudo-differential operator with the symbol — The precise definitions will be 
given below. 

In the present paper we construct the random walks the governing equation of 
which is a distributed space fractional order differential equation 

d f 2 

(2) — u(t,x) = / a(a)D%u(t, x)da, t > 0, x e R N , 



dt 

where a(a) is a positive integrable function (positively defined distribution). 

The study of properties of distributed order differential operators and their ap- 
plications to diffusion processes has been developed extensively in recent years, 
although such operators were first mentioned by Caputo @, Q in 1960th. The dis- 
tributed order differential equations have been used by Bagley, et al. Q to model 
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the input-output relationship of linear time-variant system, by Lorenzo, et al.[22j 
to study of rheological properties of composite materials, by Chechkin, et al. [8(] to 
model some ultraslow and lateral diffusion processes. Diethelm, et al. |9| studied the 
numerical aspects of such equations. Umarov, et al. [37| studied general properties 
of distributed order differential equations and solvability problems of the Cauchy 
and multipoint value problems. 

The method used in this paper for construction of multi-dimensional random 
walks are essentially based on the symbolic calculus of pseudo-differential operators 
and on the convergence properties of some simple cubature formulas. This method 
is new even in the one-dimensional case and was suggested recently in 35, [36j]. 
We note that the scaling limits are obtained in terms of characteristic functions 
of transition probabilities. The equivalence of corresponding convergence notions is 
well-known (see, [THEl). See also the recent book by M. Meerschaert and SchcfHer 
where the multi-dimensional operator stable probability distributions are studied 
and analogs of different type limits considered. Multi-dimensional random walks are 
frequently used in modeling various processes in different areas PLUM mm US- 

The present report is organized as follows. In Section 2 we give preliminaries 
simultaneously introducing the terminology that will be used in the paper. We also 
recall some properties of pseudo-differential operators with constant symbols and 
lay out some elementary properties of symbols. These properties play an essential 
role later in the study of the diffusion limits of random walks. In Section 3 we 
formulate our random walk problem in terms of sequences of i.i.d. (independent 
identically distributed) random vectors. In this Section we also formulate the main 
results. 



2. Preliminaries 

We use the following notation. is the iV-dimensional Euclidean space with 
coordinates x = (xi, . . . , xn)', Z n is the TV-dimensional integer- valued lattice with 
nodes j = ( ji , . . . , jn ) ■ We denote by xj = [hj\ , . . . , hjjy), j G 7L N , the nodes of the 
uniform lattice Zj^ defined as (KE}j N with a positive number h, the mesh width. 

We assume that a walker is located at the origin xq — (0, . . . , 0) at the initial time 
t = 0. In our random walk, at every time instants t\ = T,t% = 2r, . . . , t n — tit, . . . 
the walker jumps through the nodes of the lattice Z^. By pj, j 6 Z N , we denote 
transition probabilities. Namely, pj means a probability of jumping of the walker 
from a point Xk £ to a point Xj+k £ Zj^, where j and k are in Z w . Transition 
probabilities satisfy the non-negativity and normalization conditions: 

1- Pj>0,jeZ N ; 

Transition probabilities {pj,j G Z N } are associated with a discrete function p : 
7L N [0, 1 } . For given two transition probabilities, p and q we define the convolution 
operation p * q by 

fcez™ 

Let / be a continuous function integrable over R^. Then, as is known the 
rectangular cubature formula 



(3) 
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is valid. 

The operators in our random walk models have a close relationship to pseudo- 
differential operators with the symbols depending only on the dual variables. Sym- 
bols are allowed to have singularities. For general orientation to the theory of such 
operators we refer to [10, 0, lla- UM. l20l. l2ll . 

Let A(D), D = (£)i, . . . , D^), Dj = d/idxj,j = 1, . . . , TV, be a pseudo-differential 
operator with a symbol A(£) not depending on x, and defined in l w . We refer to 
the variable £ as a dual variable. Both type of variables, x and £ belong to R (more 
precisely, £ belongs to the conjugate (1^)" = M N ). To avoid confusion sometimes 
we write M.^ and R^, indicating their relevance to the variables x and £ respec- 
tively. Further, for a test function f(x) taken from the classical space S(R^) of 
rapidly decreasing functions, the Fourier transform 

is well defined and belongs again to S(M.¥). Let S'(M. N ) be the space of tempered 
distributions, i.e. the dual space to S(M. ). The Fourier transform for distributions 
/ 6 S"(R^) is usually defined by the extension formula (/(£), vKO) = <fi( x ))i 
with the duality pairing (., .) of S'(Rf ) and 5(Rf ). 

Definition 2.1. Assume G to be an open domain in R^. Let a function / be 
continuous and bounded on R^ and have a Fourier transform (taken in the sense 
of distributions) /(£) with compact support in G. We denote by \E'g(R^) the set of 
all such functions endowed with the following convergence. A sequence of functions 
fm £ ^g(R^) is said to converge to an element fo G \E , G(R a; Ar ) iff: 

1. there exists a compact set K C G such that supp f rn C K for all m — 1,2,...; 
2- ||/m - /oil = sup \f m - f \ — > for to -> oo. 

In the case G = R^ we write simply ^(R X N ) omitting R^ in the index of *g(R^)- 

Note that according to the Paley- Wiener theorem functions in ^g(R^) are entire 
functions of finite exponential type (see [53|, (lol|). 

Denote by H S (R% ),s G (-oo,+oo) the Sobolev space of elements / G S"(Rf ) 
for which (1 + |£| 2 ) s/2 |/(£)| G L 2 (Rf )■ It is known that if / G L„(K/) with 
p > 2, then its Fourier transform / belongs to H~ S (M.^), s > iV(| — ~). Letting 

p -> oo we get / G iJ" s (Rf ),s > f for / G L oc (R r JV ). It follows from this fact 
and the Paley- Wiener theorem that the Fourier transform of / G \E , g(R jv ) belongs 
to the space 

where H~ S (G) is a negative order Sobolev space of functionals with compact sup- 
port on G. Hence / is a distribution, which is well defined on continuous functions. 

Let ^_ G (R N ) be the space of all linear bounded functionals defined on the space 
\E , g(R JV ) endowed with the weak (dual with respect to , J , g(R JV )^) topology. By the 
weak topology we mean that a sequence of functionals g m G *S?_ G (M. N ) converges 
to an element go G v I'_ G (R Ar ) in the weak sense if for all / € ^G(R Ar ) the sequence 
of numbers (g m , /} converges to (go, f) as to — > oo. By (g, f) we mean the value of 
g G V'_ G (R N ) on an element / G * G (R JV ). 
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Definition 2.2. Let A(t;) be a continuous function denned in G C MP. A pseudo- 
differential operator A(D) with the symbol A(£) is defined by the formula 

(4) A(DMx) = (0, A(0e- ix % V G *a(R N ). 

Obviously, the function A{^)e~ lx ^ is continuous in G. Thus, A(D) in Eq. (H]) is 
well defined on f ^(R^). If <p is an integrable function with supp<p C G, then ([4]) 
takes the usual form of pseudo-differential operator 

A(D)cp(x) = J A(0me-^d4, 

with the integral taken over G. Note that in general this integral may not make 
sense even for infinitely differentiable functions with finite support (see |14|). 

Now we define the operator A(—D) acting in the space *ff_ G (WL N ) by the exten- 
sion formula 

(5) < A(-D)f, <p> = <f, A(D) V >, / e *'_ G (E*), v G * G (R JV ). 

We recall (see (HI) some basic properties of pseudo-differential operators introduced 
above. 

Lemma 2.3. The pseudo-differential operators A{D) and A{—D) with a continuous 
symbol A(£) act as 

1. A{D) : V G (m. N ) -> V G (R N ), 

2. A(-£>) : *'_ G (R JV ) -> V'_ G {R N ) 

respectively, and are continuous. 

Lemma 2.4. Let be a function continuous on H N . Then for £ £ M. N 

A{D){e- txi } = A{£)e- ix(L . 
Proof. For any fixed £ £ l w the function e"" 4 is in *(R Ar ). We have 

A( J D){e—«} = ? -l 57 / A(r?)e- fa "^(»?), 

where dMf(»?) = ^[e - ^]^ = (2tt) n 6(t) - £) d ?7- Hence A{D){e-^} = A(£)e" M «. 

□ 

Corollary 2.5. 1. A(£) = (A(D) e - ra «)e'^; 

2. A(0 = (A(Z>)e--€)| ic=0 ; 

3. A(£) =< A(—D)5(x),e~ lx £ >, where S is the Dirac distribution. 

Remark 2.6. Since the function e" 43 ^ does not belong to S(R N ) and D(R N ), the 
representations for the symbol obtained in Lemma 12.41 and Corollary 12.51 are not 
directly applicable in these spaces. 

Denote by Dq , < a < 2, the pseudo-differential operator with the symbol 
— |£| a . It is evident that Dq coincides with the Laplace operator A for a — 2. For 
a < 2 it can be considered as a fractional power of the Laplace operator, namely 
Dn = —(—A)"/ 2 . Dq can also be represented as a hypersingular integral (see, e.g. 

M) 

(6) D%f(x) = 6(a) / jj^dy, 

Jr" \y\ 
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where A y is the second order centered finite difference in the y direction, and b(a) 
is norming constant defined as 

Q r(f)r(^)sin^f 

( ? ) b ( a > = 2 2- a7r l+JV/2 ■ 

It is seen from ([7]) that in the representation ([6]) of Dq the value a = 2 is singular. 
Lemma 2.7. for i/ie symbol of Dq the following equalities hold true: 

A 2 \ ixt, 



(Dge^) | x=0 - 6(a) / 1 -fedw|x=o = -|£| Q , < a < 2. 

jk« \y\ 



(8) 



Proof. This statement is a direct implication of Corollary |2 . 5l applied to the operator 
DJ. □ 

The cubature formula ([3]) yields for the integral in the right hand side of ([6]) 

L \y\ N+a y k h:»w N+a ' 

where fj — f(xj) and fc| is Euclidean norm of k — (fci, . . . , kpj) G 
Consider the distributed space fractional order differential equation 

d f 2 

(10) —u{t,x) = a{a)D o Y u{t,x)da,t>0,x£R N , 
ot J 

where a(a) is a positive (in general, generalized) function defined in (0, 2]. A distri- 
bution G(t, x), which satisfies the equation (|10p in the weak sense and the condition 

(11) G(0,x) = 5(x),x 6 R N , 

where 5(x) is the Dirac's distribution, is called a fundamental solution of the Cauchy 
problem (JTDJ), (fTTj) . In the particular case of 

M 

(12) a(a) = a m 6(a - a m ), < a x < ■ ■ ■ < a M < 2, 

m— 1 

with positive constants a m we get a multiterm space fractional differential equation 

M 

,/l / r i : 



9 

(13) — x) = 2_, a m D* m u(t, x)t> 0, x eR N . 



Denote the operator on the right hand side of the equation (JTUJ) by B(D). It can be 
represented as a pseudo-differential operator with the symbol 

(14) B(0 = - / 2 a{a)\i\ a da. 

Jo 

It is not hard to verify that the fundamental solution of equation (|10|l is 

(15) G(t, x) = F^ 1 (e tB ^ 
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where F 1 stands for the inverse Fourier transform. In the particular case of a(a) — 
8 (a — 2) we have the classical heat conduction equation 

—u(t,x) = Au(t,x), t > 0, x e R*, 

whose fundamental solution is the Gaussian probability density function evolving 
in time 

For a{a) — 8{a — ao), < ao < 2, the corresponding fundamental solution is the 
Levy ao-stable probability density function [301 ] 

(16) Gaoi t,x) = -± w [ e-'l«l Q V*«de 

(27T) JV J r n 

The power series representation of the stable Levy probability density functions is 
studied in [3J, [23J, |33j . Recall also that ao = 1 corresponds to the Cauchy-Poisson 
probability density (see [28j]) 

Gx(t,x)- [ 2 ' 



7r (n+l)/2 (j x \2 +i 2)(n+l)/2 



3. Main results: construction of random walks 



In this Section we construct random walks associated with distributed space frac- 
tional order differential equations (| 10[) . More precisely, we show that the special 
scaling limit of the constructed random walk is a diffusion process whose probabil- 
ity density function is the fundamental solution of (|10j) . 

Let X be an N-dimensional random vector [25| which takes values in Z N . Let 
the random vectors Xi , X2 , . . . also be N-dimensional independent identically dis- 
tributed random vectors, all having the same probability distribution, common with 
X. We introduce a spatial grid {xj = jh,j € Z N }, with h > and temporal grid 
{t n = nr, n = 0,1,2,...} with a step r > 0. Consider the sequence of random 
vectors 

S n = hX-i + /1X2 + • • • + /iX„, n = 1, 2, . . . 

taking So = for convenience. We interpret Xi, X2, . . . , as the jumps of a particle 
sitting in x = xq = at the starting time t = to = and making a jump X n from 
S„_i to S„ at the time instant t = t n . Then the position S(t) of the particle at 
time t is 

l<k<t/r 

Denote by Uj(t n ) the probability of sojourn of the walker at Xj at the time t n . 
Taking into account the recursion S n _j_x — S n + /iX„ we have 

(17) Vj{tn+l)= E PkVj-k(t n ),j ^ Z N ,Tl = 0, 1, . . . 

fcez™ 

The convergence of the sequence S„ when n — > oo means convergence of the discrete 
probability law (yj(t n ))jez N ? properly rescaled as explained below, to the probabil- 
ity law with a density u(t, x) in the sense of distributions (in law). This is equivalent 
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to the locally uniform convergence of the corresponding characteristic functions (see 
for details [25J| ) . We use this idea to prove the convergence of the sequence of char- 
acteristic functions of the constructed random walks to the fundamental solution 
of distributed order diffusion equations. 

In order to construct a random walk relevant to (|10p we use the approximation 
(J3]) for the integral on the right hand side of (J6j> , namely 

fcez" 1 1 

and the first order difference ratio 

du ^ Uj(t n+1 ) - Uj(t n ) 
dt~ t 

for with the time step r = t/n. Then from PU|) we derive the relation (fTT|) with 
the transition probabilities 



(18) Pk 
where 



i-2T£ m#0 ^h^ = 0; 

2r%^, iffc^O, 



(19) Q m {h) - ^ \m\- a P (a, h)da, p(a, h) = <°^} a \ 
Assume that the condition 

(20) ^M^arEW<i. 

is fulfilled. Then, obviously, the transition probabilities satisfy the properties: 

i- Sfcez« Pk = i; 

2. p k > 0,k e Z N . 
Introduce the function 

K(a) = y = T -¥rh, < a < 2, 

v ' / ' U. Af+Q / t m N+a ' — ' 

k^O 1 1 m=l 

where M m = X)|fc|=m ^ ne one-dimensional case 7^.(a) coincides with the Rie- 

mann's zeta-function, IZ(a) — 2£(1 + a).) The Eq. (|2"D|) can be rewritten as 

(21) j(T >) = 2 r / 2 °»M ^<i. 

It follows from the latter inequality that h — > yelds t — > 0. This, in turn, yields 
rt = t/r — > oo for any finite i. 

Now we assume that the singular support of a does not contain 2, i.e., {2} ^ 
singsupp a 0- 



1 This condition relates only to distributions. 
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Theorem 3.1. Let X be a random vector with the transition probabilities pk — 
P(X = .Tfe), k 6 1 N , defined in Eqs: (j 1 8[) . (|19[) and, which satisfy the condition (|20p 
(or, the same, (|21| J. Then the sequence of random vectors S n — /iXi + • • • + ftX n , 
converges as n — > oo in law to the random vector whose probability density func- 
tion is the fundamental solution of the distributed space fractional order differential 
equation (QUI) . 

Proof. We have to show that the sequence of random vectors S„ tends to the 
random vector with pdf G(t,x) in Eq. (|15|) . or the same, the discrete function 
yj(t n ) tends to G(t,x) as n — > oo. It is obvious that the Fourier transform of 
G(t,x) with respect to the variable x is the function G(i, £) = e tB ^\ where £>(£) is 
defined in Eq. (fl4|) . Let j3(— £) be the characteristic function corresponding to the 
discrete function pj~, k G 1> N , that is 



kez, N 



It follows from the recursion formula (jTTJ) (which exhibits the convolution) and 
the well known fact that convolution goes over in multiplication by the Fourier 
transform, the characteristic function of yj{t n ) can be represented in the form 

Taking this into account it suffices to show that 

(22) P n {-K) -> e tB ^\n -> oo. 

The next step of the proof is based on the following simple fact: if a sequence s n 
converges to s for n — > oo, then 

(23) lim(l + — ) n = e s . 

n 

We have 

p«(-^) = (l-r^^(l-e^)r 
Mo 11 



^ i*r y 



It follows from ([3]) and Corollary 12 . 51 that 



/\2 e ik£h 

^ \kh\ N + a ' 



& ( tt ) V , V v 



tends to (Dq e tx £)i = — as ft — > (or, the same, n —* oo ) for all a 6 (0,2] 
Hence 

s n = J a(a){b(a) £ __-ft^}cfo -> n - oo (ft - 0). 



Thus, in accordance with ([23]) we have 

p"(-ftO e w ^\ n -> oo. □ 
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The random walk related to the multiterm fractional diffusion equation can be 
derived from Theorem 13.11 Assume that a(a) has the form (p~2|) with < a% < 
■ ■ ■ < ccm < 2. So, we again exclude the case {2} £ singsupp a. 

Theorem 3.2. Let the transition probabilities pk = P(X = Xk),k 6 Z w , of the 
random vector X be given as follows: 



Pk 



j^O 11 m=l |J| 



1 m=l ' 



z/fc^O, 



where fi m = 7^77 m = 1, . . . , M. Assume, 



A I 



^ a m b(a m )TZ(a m )fi m < 1. 



T/ien </ie sequence of random vectors S„ = ft,Xi + • • - + /iX n , converges as n — > oo 
m Zow to t/ie random vector whose probability density function is the fundamental 
solution of the multiterm fractional order differential equation h!3\l . 



Remark 3.3. The condition {2} ^ singsupp a is required due to singularity of the 
value a = 2 in the definition of Dq (see ([?]))■ The particular case a(a) = S(a — 2) 
reduces to the classic heat conduction equation and corresponding random walk is 
the classic Brownian motion. In more general case of a(a) = Y^iiLo ci8^(a — 2) this 
condition leads to the scaling limit with cr(r, h) — h 2 lnj^ (see, also [16|). 

This work was supported in part by NIH grant P20 GM067594. 
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